%% Note: this is a file called by Main_iid_calibration or whenever it needs recalibration

%%%% We will calibrate the model (some parameters are given).
%%%% The calibration go through steps below.

%% Step 1: calibrate CHI_PI, CHI_q, XI, and G; ALPHA_0 will be used later
disp('--------------> Step 1 calibration <--------------')
disp('I am calibrating the basic parameters...')

cali_step   = 1; %%%% this indicates the calibration does not take into account acquisition elasticity

x0 = [par.EPS_L_BOUND * 4 + 2, 0.08, 0.1, par.CHI_PI, par.CHI_q, par.XI];
x  = fsolve(@(x) eqm_ss_iid_case_6(x, par, fun, 1, cali_step), x0);

[fval, eqm] = eqm_ss_iid_case_6(x, par, fun, 1, cali_step);

disp('Values of L, B, ALPHA_0, CHI_PI or CHI_k, CHI_q, CHI_XI')
disp(x);
L                   = x(1);
B                   = x(2);
par.ALPHA_0         = x(3);
switch par.credit_parameter
    case 'CHI_PI'
        par.CHI_PI  = x(4);
    case 'CHI_k'
        par.CHI_k   = x(4);
end
par.CHI_q           = x(5);
par.XI              = x(6);
par.G               = eqm.G;
par.CHI_PI_k        = par.CHI_PI;
par.surplus_net     = eqm.surplus_net;
%- par.iota * eqm.Z / eqm.w * par.XI;


par.aqc_spending    = eqm.E_a_K  / par.ALPHA_0 / eqm.prob_A;

disp('Investment-output ratio (note: this calibrates par.ETA)')
disp((par.DELTA * eqm.K )/ eqm.Y * (1 - par.IM_share))


%% Step 2: calibrate the entry-cost distribution by using ALPHA_0
disp('--------------> Step 2 calibration <--------------')

cali_step               = 2; %%%% now the calibration takes into account acquisition elasticity
par_sensitivity         = par;
par_sensitivity.iota    = par.iota + par.eps_;

if par.THETA == 0.50
    disp('I am calibrating the entry cost parameters...')
    c_mu        = par.c_mu;
    x0          = [L, B, par.ALPHA_0, c_mu]; % initial guess
    [x, ~]      = fsolve(@(x) eqm_ss_iid_case_6(x, par_sensitivity, fun, 1, cali_step), x0);
    c_mu        = x(4);
    c_sig       = (log(par.surplus_net) - c_mu) / erfinv(2 * par.ALPHA_0 - 1) / sqrt(2); 

else
    % minimize distance to the target elasticity;
    d           = 0.5; % arbitrary; but it turns out not matter much
    c_mu_ini    = -1.5;
    x0          = [L, B, par.ALPHA_0, c_mu_ini]; % initial guess        
    c_mu_grid   = [c_mu_ini + abs(c_mu_ini) * d, log(par.surplus_net) + 0.002];
    elasticity  = zeros(1, length(c_mu_grid));
    
    for ii_mu   = 1:length(c_mu_grid)
        c_mu      = c_mu_grid(ii_mu);
        c_sig     = (log(par.surplus_net) - c_mu) / erfinv(2 * par.ALPHA_0 - 1) / sqrt(2);
        par_sensitivity.c_mu  = c_mu; 
        par_sensitivity.c_sig = c_sig;
        x00       = [L, B, par.ALPHA_0];
        x         = fsolve(@(x) eqm_ss_iid_case_6(x, par_sensitivity, fun, 1, 0), x00);
        [fval, eqm] = eqm_ss_iid_case_6(x,  par_sensitivity, fun, 1, 0);
        ALPHA_0 = x(3);
        elasticity(ii_mu) = (eqm.E_a_K  /  (ALPHA_0 * eqm.prob_A) - par.aqc_spending) / par.aqc_spending / par_sensitivity.iota_per_change;
    end

    if elasticity(1) - par.aqc_elasticity_target < 0 % so that you can hit the elasticity target
        [x,fval]  = fsolve(@(x) eqm_ss_iid_case_6(x, par_sensitivity, fun, 1, cali_step), x0);
        c_mu      = x(4);
    else
        c_mu = c_mu_grid(1);
    end
    c_sig     = (log(par.surplus_net) - c_mu) / erfinv(2 * par.ALPHA_0 - 1) / sqrt(2); % at the calibrated par.ALPHA_0 and par.surplus_net
end


%% Step 3:
fprintf('--------------> Step 3: Summary <-------------------- \n')
eqm.IM_s      = average_log_normal(par.surplus_net, c_mu, c_sig);
eqm.IM        = eqm.IM_s / par.XI;


%%%% wage-pemium approach (XI_IM = XI, but the wage is different for financial services)
par.XI_IM = par.XI;
par.ZETA  = par.IM_share * eqm.Y / eqm.IM / eqm.w; % wage premium parameter discussed in the main text
par.c_mu  = c_mu; 
par.c_sig = c_sig; 

%%%% different XI approach (XI_IM and XI are different, but the wage rate
%%%% is the same for financial services (not used)   
% par.ZETA = 1;
% par.XI_IM = w * IM_s  / (par.IM_share * Y);

par.W_calibrated    = log(eqm.c_star) - par.XI * eqm.H - par.XI_IM * eqm.IM;
fprintf('the credit parameters are (chi_pi, chi_k, and chi_q): \n %.4f, %.4f, %.4f \n\n',...
    par.CHI_PI, par.CHI_k, par.CHI_q )
fprintf('The mean and the std.dev. of (log) cost of participation \n %.4f, %.4f \n\n',par.c_mu, par.c_sig)
